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We investigate the computational requirements for all-sky, all-frequency searches for gravitational 
waves from spinning neutron stars, using archived data from interferometric gravitational wave de- 
tectors such as LIGO. These sources are expected to be weak, so the optimal strategy involves coherent 
accumulaton of signal-to- noise using Fourier transforms of long stretches of data (months to years). 
Earth-motion-induced Doppler shifts, and intrinsic pulsar spindown, will reduce the narrow-band 
signal-to-noise by spreading power across many frequency bins; therefore, it is necessary to correct 
for these effects before performing the Fourier transform. The corrections can be implemented by 
a parametrized model, in which one does a search over a discrete set of parameter values (points 
in the parameter space of corrections). We define a metric on this parameter space, which can be 
used to determine the optimal spacing between points in a search; the metric is used to compute 
the number of independent parameter-space points N p that must be searched, as a function of ob- 
servation time T. This method accounts automatically for correlations between the spindown and 

P-h ' Doppler corrections. The number N P (T) depends on the maximum gravitational wave frequency and 

the minimum spindown age r = /// that the search can detect. The signal-to-noise ratio required, 

/-o ' in order to have 99% confidence of a detection, also depends on N P (T). We find that for an all-sky, 

all-frequency search lasting T = 10 7 s, this detection threshhold is h c fa (4 5)h s / yr , where h 3 / yl 

is the corresponding 99% confidence threshhold if one knows in advance the pulsar position and spin 
£> period. 

We define a coherent search, over some data stream of length T, to be one where we apply a correction, 
followed by an FFT of the data, for every independent point in the parameter space. Given realistic 
limits on computing power, and assuming that data analysis proceeds at the same rate as data 
acquisition (e.g. 10 days of data gets analysed in ~ 10 days), we can place limitations on how much 
data can be searched coherently. In an all-sky search for pulsars having gravity-wave frequencies 
/ < 200Hz and spindown ages r > lOOOYrs, one can coherently search ~ 18 days of data on a 
teraflops computer. In contrast, a teraflops computer can only perform a ~ 0.8-day coherent search 
for pulsars with frequencies / < 1kHz and spindown ages as low as 40Yrs. 

In addition to all-sky searches we consider coherent directed searches, where one knows in advance 
5-H ■ the source position but not the period. (Nearby supernova remnants and the galactic center are 

obvious places to look.) We show that for such a search, one gains a factor of ~ 10 in observation 
time over the case of an all-sky search, given a lTflops computer. 

The enormous computational burden involved in coherent searches indicates a need for alternative 
data analysis strategies. As an example we briefly discuss the implementation of a simple hierarchical 
search in the last section of the paper. Further work is required to determine the optimal approach. 
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I. INTRODUCTION 

The direct observation of gravitational waves is a re- 
alistic goal for the kilometer-scale interferometers which 
are now under construction at various sites around the 
world jy,|j. However, the battle to see these waves is not 
over when the detectors are constructed and running. 
Searching for gravitational wave signals in the interfer- 
ometer output presents its own problems, not least of 
which is the sheer volume of data involved. 

Potential sources of gravitational waves fall roughly 
into three classes: bursts, stochastic background, and 
continuous emitters. 

Burst sources produce signals which last for times con- 
siderably shorter than available observation times. The 
chirp signals from compact coalescing binaries belong to 
this class. Since theoretical waveforms, valid during the 
inspiral phase of the binary evolution, have been accu- 
rately calculated using post-Newtonian methods J3|, it 
is possible to search the data stream for chirps using 
matched filtering techniques. Detailed studies have been 
carried out to ascertain the optimal set of search tem- 
plates Qg], and a preliminary investigation of search 
algorithms is now under way ||. Detection of other, 
not so well understood, sources in this class — e.g. non- 
axisymmetric supernovae — has received limited atten- 
tion§. 

Flanagan || has determined how to cross-correllate the 
output of two detectors in order to search for a stochastic 
background of gravitational radiation, which was imple- 
mented by Compton |9) and applied to data taken during 
a period of 100 hours by two prototype interferometers 
detectors in Glasgow and Garching |l0[ . In pl[ , Allen 
presents a detailed discussion of the potential significance 
of detecting a stochastic background. Compton's work, 
and simulations performed by Allen, have demonstrated 
that this kind of analysis requires minimal computational 
resources. 

In this paper we consider some issues involved in 
searching for continuous wave sources. Throughout our 
discussion we use pulsars as a guide to develop a search 
strategy. 



A. Gravitational waves from pulsars 

Rapidly rotating neutron stars (pulsars) tend to be ax- 
isymmctric; however, they must break this symmetry in 
order to radiate gravitationally. The pulsar literature 
contains several mechanisms which may lead to deforma- 
tions of the star, or to precession of its rotation axis, and 
hence to gravitational wave emission. The characteristic 



amplitude^ of gravitational waves from pulsars scales as 

h c ~ ^ (1.1) 

r 

where / is the moment of inertia of the pulsar, / is the 
gravitational wave frequency, e is a measure of the de- 
viation from axisymmetry and r is the distance to the 
pulsar. 

Pulsars are thought to form in supernova explosions. 
The outer layers of the star crystallize as the newborn 
pulsar cools by neutrino emission. Estimates, based on 
the expected breaking strain of the crystal lattice, suggest 
that anisotropic stresses, which build up as the pulsar 
loses rotational energy, could lead to e < 10~ 5 ; the exact 
value depends on the breaking strain of the neutron star 
crust as well as the neutron star's 'geological history', and 
could be several orders of magnitude smaller. Nonethe- 
less, this upper limit makes pulsars a potentially interest- 
ing source for kilometer scale interferometers. Figure H 
shows some upper bounds on the amplitude due to these 
effects. 

Large magnetic fields trapped inside the superfluid in- 
terior of a pulsar may also induce deformations of the 
star. This mechanism has been explored recently in p2| , 
indicating that the effect is extremely small for standard 
neutron star models (e < 10 -9 ). 

Another plausible mechanism for the emission of 
gravitational radiation in very rapidly spinning stars 
is the Chandrasekhar-Friedman-Schutz (CFS) instabil- 
ity, which is driven by gravitational radiation reac- 
tion |[3uL4j. It is possible that newly- formed neutron 
stars may go through this instability spontaneously as 
they cool soon after formation. The radiation is emit- 
ted at a frequency determined by the frequency of the 
unstable normal mode, which may be less than the spin 
frequency. 

Accretion is another way to excite neutron stars into 
emitting gravitational waves. Wagoner [ |15[ proposed 
that accretion may drive the CFS instability. There is 
also the Zimmermann-Szcdinits mechanism mM where 
the principal axes of the moment of inertia are driven 
away from the rotational axes by accretion from a com- 
panion star. Accretion can in principle produce relatively 
strong radiation, since the amplitude is related to the ac- 
cretion rate rather than to structural effects in the star. 
However, accreting neutron stars will be in binary sys- 
tems, and these present problems for detection that go 
beyond the ones we discuss in this paper. We hope to re- 
turn to the problem of looking for radiation from orbiting 
neutron stars in a future publication. 



1 We adopt the definition of h c provided in Eq. (50) of 
Thome iffll. 



B. Three classes of sources 

Observed pulsars fall roughly into two groups: (i) 
young, isolated pulsars having periods of tens or hun- 
dreds of milliseconds, and (ii) older, millisecond pulsars. 
The young pulsars are most likely to deviate significantly 
from axisymmetry; however, they are generally observed 
to have low frequencies, so that there is a competition be- 
tween the frequency, /, and deviation from axisymmetry, 
e, in Eq. (1.1). On the other hand, millisecond pulsars, 
whose waves are higher in frequency, tend to be quite old 
and well annealed into an axisymmetric configuration. 

Radio observations can only probe a small portion of 
our galaxy in searching for pulsars. A significant ef- 
fect reducing the depth of radio searches is dispersion of 
the signal by galactic matter between potential sources 
and the earth. Given current evolutionary scenarios for 
pulsars — that they are born in supernova explosions — 
it seems likely that most pulsars should be located in 
the galactic disk, and the youngest of these will also be 
shrouded in a supernova remnant, making them invisible 
to radio astronomers. 

Blandford |17|,f7| has pointed out that there could ex- 
ist a class of pulsars which spin down primarily due to 
gravitational radiation reaction. For sources in this class 
the frequency scales as / ex r -1 / 4 , where r is the age 
of the pulsar. If the mean birth rate for such pulsars in 



It is important that any search strategy should be gen- 
eral enough to encompass all three of the above classes, 
allowing for the significant changes in frequency which 
may be inherent in the sources (see section ||) . 
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FIG. 1. Characteristic amplitudes h c [see Eq. (3.5)] for 
several postulated periodic sources, compared with sensitivi- 



, • — 1 .I i it 1 ]• . ties /i-i/ vr of the initial and advanced detectors in LIGO. (h-M 

our galaxy is r R , the nearest one should be a distance ^/yr ,.,,,, , , 



r = R^tb/t from earth, where R ~ lOkpc is the radius 
of the galaxy. The intrinsic gravitational wave ampli- 
tude (that is, the amplitude h at some fixed distance) of 
a pulsar in this class is proportional to t~ 1//2 . Thus, the 
nearest source in this class would have a dimcnsionlcss 
amplitude h c at the Earth 
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(1.2) 



independent of the frequency, or the the ellipticity e of 
the pulsar. Assuming the existence of such a class of 
pulsars, with tb < 2 x 10 4 Yrs, we see from Fig. [l] that 
there is a large region of parameter space that is both 
(i) detectable by the LIGO detector and (ii) physically 
reasonable, in the sense that e < 10~ 5 and / lies in the 
range 200-1000 Hz. 

Note that Blandford's argument can be slightly re-cast 
to yield an upper limit on the gravitational wave strength 
of any isolated pulsar-i.e., any pulsar whose radiated an- 
gular momentum is not being replenished by accretion. 
The age of an isolated pulsar must be shorter than the 
age computed assuming the spin-down is solely due to 
gravitational wave emission. Correspondingly, if we set 
t b equal to 40 yrs (corresponding to the birthrate for 
all pulsars), we get the following upper limit for mea- 
sured gravitational wave amplitude of an isolated pulsar: 
h c < 2 x 10~ 24 . Of course, this is a statistical argument. 
This bound could certainly be violated by an isolated 
pulsar that just happens to be anomalously close to us. 



corresponds to the amplitude h c of the weakest source de- 
tectable with 99% confidence in |yr = 10 7 s integration time, 
if the frequency and phase of the signal is known in advance.) 
Long-dashed lines show the expected signal strength as a func- 
tion of frequency for pulsars at a distance of 10 kpc, assuming 
non-axisymmetries of e = 10~ 5 and e = ICE 8 , where e is de- 
fined in section |TJA. Upper limits are also plotted for the 
Crab and Vela pulsars, assuming their entire measured spin- 
down is due to gravitational wave emission. The dotted lines 
indicate the strongest waves received at the earth for Bland- 
fords hypothetical class of pulsars; each line corresponds to a 
particular birth rate. 



C. The data analysis problem 

The detection of continuous, nearly fixed frequency 
waves will be achieved by constructing power spec- 
trum estimators and searching for statistically significant 
peaks at fixed frequencies. In practice, this is achieved by 
calculating the amplitude of the Fourier transform of the 
detector output given by applying a fast Fourier trans- 
form (FFT) , a discrete approximation to the true Fourier 
transform: 



Hf) 



1 



e 2 " /4 h(t) dt . 



(1.3) 



The main hope of detection lies in the fact that one may 
observe the sky for long time periods of time T. When 
such a data stretch is transformed to make the underlying 



signal monochromatic, the signal to noise ratio grows as 
VT in amplitude (or as T in the power spectrum). One 
will likely need to have integration times of several weeks 
or months in order for the expected signals from nearby 
sources to rise above the noise. However, such long data 
stretches pose a significant computational burden; using 
10 7 s of data to look for signals with gravitational wave 
frequencies up to 500Hz requires calculating an FFT with 
N ~ 10 10 data points. Calculation of a single such FFT 
would take about Is on a lTflops computer, assuming 
that all 10 10 points can be held in fast memory. Unfor- 
tunately, this is not the whole story. 

The detection problem is complicated by the fact 
that the signal received at the detector is not perfectly 
monochromatic. Earth-bound detectors participate in 
complex motions which lead to significant Doppler shifts 
in frequency as the Earth rotates, and as it orbits around 
the sun (this orbit is significantly perturbed by the moon 
and the other planets). The time dependent accelera- 
tions broaden the spectral lines of fixed frequency sources 
spreading power into many Fourier bins about the ob- 
served frequency. In order to maintain the benefit of 
long observation times, it is therefore necessary to remove 
the effects of the detector motion from the data stream. 
This can achieved by introducing an inertial (barycen- 
tered) time coordinate and carrying out the FFT with 
respect to it. The difficulty of doing this was first esti- 
mated by one of us |Ug] . However, we must also consider 
the additional complication that the signal may not be 
intrinsically monochromatic. If the signal exhibits intrin- 
sic frequency drift, or modulation, due to the nature and 
location of the source — as is expected for pulsars which 
spin down with time — these effects can also be removed 
in the transformation to the new time coordinate. 

Unfortunately, the demodulated time coordinate de- 
pends strongly on the direction from which the signal is 
expected, and on the intrinsic frequency evolution one 
assumes for the source. Thus, in searching for sources 
whose position and timing are not well known in advance 
one must apply many different corrections to the data, 
performing a new FFT after each correction. Given the 
possibility that the strongest sources of continuous grav- 
itational waves may be electromagnetically invisible or 
previously undiscovered, an all sky, all frequency search 
for such unknown sources is of considerable interest. To 
obtain some idea of the magnitude of this task, consider 
searching the entire sky for signals with (fixed) frequen- 
cies up to 500Hz using 10 7 s worth of data. Assuming the 
entire data stream could be held in fast memory on a ma- 
chine capable of lTflops, it would take 10 8 s to complete 
the search. Introducing intrinsic spindown effects into 
the search increases the computational cost, at fixed in- 
tegration time, by many orders of magnitude. This com- 
putational cost is the central problem of searching for 
unknown pulsars in the output from gravitational wave 
detectors and is the focus of this paper. 



D. Summary of results 

We parametrise the space of pulsar signals by the po- 
sition of the source on the sky {9, <fi}, entering through 
Doppler shifts due to the detector's motion, and by spin- 
down parameters /& which cha ract erise the intrinsic fre- 
quency evolution. [See Eq. ( |3.7| ).1 We constrain the 
range of possible values of the spindown parameters us- 
ing the (spindown) age r = /// of the youngest search 
that a search can detect, thus \fk\ < r~ k . For the 
computationally-intensive search over all sky positions 
and spindown parameters, it is important to be able to 
calculate the smallest number of independent parameter 
values which must be sampled in order to cover the entire 
space of signals. We have accomplished this by introduc- 
ing a distance measure and corresponding metric on the 
parameter space. The analysis is patterned after a sim- 
ilar one developed by Owen H for gravitational waves 
from inspiralling, compact binaries. Using our metric 
one can compute the volume of parameter space, thus 
inferring the number of independent points that must be 
sampled in order to cover the entire space. We define a 
coherent search to be one where we perform one demod- 
ulation and FFT of the data for every independent point 
in the parameter space. Besides telling us the compu- 
tational requirements for a coherent search, the metric 
approach tells us how to place the points most efficiently 
in parameter space, in a similar way to that discussed by 
Owen. 

We have found it useful to present the results based on 
several possible search strategies, which cover different 
regions of the parameter space. Accordingly, we define a 
pulsar to be old if its spindown age r is greater than I0 3 
Yrs and young if r > 40 Yrs. A pulsar is considered to 
be slow if its gravitational wave frequency is / < 200 Hz 
and fast if/ < 10 3 Hz. 

A coherent all-sky search of 10 7 seconds of data for old, 
slow pulsars requires approximately 1.1 x 10 10 indepen- 
dent points in the parameter space; only one spindown 
parameter is needed to account for intrinsic frequency 
evolution. In contrast, an all-sky search for fast, young 
pulsars in 10 7 seconds of data requires 8 x 10 21 indepen- 
dent parameter space points to be sampled, using three 
spindown parameters to model intrinsic frequency evo- 
lution. Note that searches for old, fast pulsars (such as 
known millisecond radio pulsars) and young, slow pul- 
sars (younger brothers of the Crab and Vela) are au- 
tomatically subsumed under the latter search. These 
results mean the following. Assuming unlimited com- 
puter power and stationary, gaussian statistics, a pul- 
sar with unknown position and period must have strain 
h c rs 4.3/i3/ yr , if it is in our 'old, slow' category, and 
h c ~ 5.1/i 3 / yr , if it is in our 'young, fast' category, to 
be detected with 99% confidence in a 10 7 second search. 
Here h 3 / yr is the strain required for detection with 99%- 
confidcncc in a 10 second integration, assuming the pul- 



sar position and period are known in advance 



*3/yr 



(/) = 4.2 y/S n (f) x 10-^Hz 



(1.4) 



Thus, when considering an all-sky, all-frequency pulsar 
search, the LIGO sensitivity curves shown in Fig.l effec- 
tively overestimate the detector's sensitivity by a factor 
of ~ 4 5, even in the limit of infinite computing power. 

Our ability to perform searches for continuous waves 
will certainly be limited by the available computing re- 
sources. Assuming realistic computer power — say of 
order 10 13 flops — we estimate that computing limita- 
tions will effectively reduce the sensitivity of the detector 
by another factor of ~ 2, even for some reasonably opti- 
mized and efficient search strategy. However more work 
will be needed to develop a reasonably optimized algo- 
rithm, and thus to refine this latter estimate. 

While the concept of the metric is introduced in the 
framework of an all-sky search for unknown pulsars, it is 
clear that we may use the same approach to examine the 
depth of a search over limited regions of the parameter 
space. In particular, once the scope of a search is decided, 
the optimization procedure dicussed in section VI can be 
used to determine the observation time and grid spac- 
ing which maximises the expected sensitivity of a search. 
As an example, we consider coherent directed searches, 
in which one assumes a specific sky position (such as 
a particular cluster or supernova remnant) and searches 
only over spindown parameters. Again, we present re- 
sults for two concrete scenarios based on fast, young pul- 
sars and old, slow pulsars. Similar considerations apply 
to directed searches as to all-sky searches; that is, the 
curves in Fig. |l] overestimate the detector sensitivity for 
10 7 second integration. Table 1 summarises the results 
for both cases. 

We note that in each type of search, the number of 
parameter space points, and hence the computational re- 
quirements, were reduced significantly by the assumption 
that the points were placed with optimal spacings given 
by the metric formalism. Nevertheless, the bottom line is 
that limitations on computational resources will severely 
restrict the integration times that can be achieved. As- 
suming access to a Tflops of computing power (effective 
computational throughput, ignoring possible overheads 
due to interprocessor communication or data access), we 
find the following limits on coherent integration times: 
For young, fast pulsars we are limited to about 0.8 days 
for an all-sky search, and 18 days for a directed search. 
For older, slower pulsars, on the other hand, we are only 
limited to 9 days for an all-sky search, and nearly 160 
days for a directed search. The threshold sensitivities 
that these strategies can achieve, relative to the noise 



curves in Fig 
power in Fig. 



are plotted as functions of computing 



TABLE I. The number of independent parameter points 
N P (T, /imax = 0.3) required for a coherent T — 10 7 s search, for 
four fiducial types of pulsar. We list the requirements both 
for all-sky searches and for directed searches (i.e., searches 
where the source position is known in advance). Also listed 
are the threshhold values /i t h of the characteristic strain h c 
required to have 99% confidence of detection, assuming un- 
limited computer power . These thresho ld values are given by 
hth/h 3/yr = (1/1.90) V / I r i(5bA 7 7vpj -1 where N = 2/ max T. 
Here h 3 / yr is the corresponding threshhold, assuming the pul- 
sar's postion and period and are known in advance. 



Search Parameters N p 
f (Hz) r (Yrs) (all-sky) 



ftth/ft 



3/yr 



N„ 



hth/h 



3/yr 



< 200 

< 10 3 

< 200 

< 10 3 



> 10 3 

> 10 3 
>40 
>40 



(all-sky) (directed) (directed) 



1.1 x 10 1 

1.3 x 10 1 

1.7 x 10 1 

8 x 10 21 



3.7 
4.2 
4.3 
4.6 



3.7 x 10" 
1.2 x 10 8 
8.5 x 10 l; 
1.4 x 10 11 



3.3 
3.5 
3.9 

4.1 



2 This differs from equation (112) in |7j because we have spec- 
ified 99% confindence, and we have use the correct exponential 
probability function for power. 
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FIG. 2. Relative amplitude sensitivities h, 

achievable with given computational resources 
ous coherent search strategies: (a) directed search for old 
(t > lOOOYrs), slow (/ < 200Hz) pulsars, (b) all-sky search 
for old, slow pulsars, (c) directed search for young (r > 40Yrs) 
fast (/ < 1000Hz) pulsars, and, (d) all-sky search for these 
same sources, hi/r is the expected sensitivity of the detector 
for an observation time T seconds, and with 99% confidence, 
assuming only that the frequency bandwidth of the source is 
constrained in advance. For a given computational power, we 
have determined the opt imum observation time as described 
in sections [VIB and VII . 



E. Organization of this paper 

In section |l| we outline the physics of pulsars which 
is relevant to the detection of continuous gravitational 
waves. The discussion is phenomenological and based 
almost entirely on pulsar data collected by radio as- 
tronomers. We focus attention on effects which may lead 
to significant frequency evolution over periods of several 
weeks of observation. 

Then, in section 



III, we introduce a parameterised 



model of the expected gravitational waveform, including 
modulating effects due to detector motion. 

From this, we go on in section |y| to describe the ba- 
sic technique used to search for signals, by constructing 
a demodulated time series. Livas |l9[ , Jones [E0[ and 
Niebauer [gl| have implemented variants of this basic 
search strategy over limited regions of parameter space 
(in particular they have not considered pulsar spin-down, 
and have restricted attention to small areas of the sky) . 



For the more computationally-intensive search over all 
sky positions and spindown parameters, it is important 
to be able to calculate the smallest number of indepen- 
dent parameter values which must be sampled in order 
to cover the entire space of signals. In section M we de- 
velop the metric formalism for calculating the number of 
independent points in para meter space. 

In sections VI and VII we apply this formalism to 
determine the computational requirements of an all-sky 
search for unknown pulsars and a directed search, respec- 
tively. 



Finally in section VIII, we list some possible alterna- 
tives to a straightforward coherent search of the inter- 
ferometer data. Detailed studies of the pros and cons of 
each are currently under investigation. 



II. PULSAR PHENOMENOLOGY 

Currently, the only expected sources of continuous, pe- 
riodic gravitational waves in the LIGO band are pulsars. 
In this section, therefore, we review those properties of 
pulsars which may be important in the detection pro- 
cess. In general, the search technique we present later 
is capable of detecting any nearly monochromatic grav- 
itational wave with sufficient amplitude. However, it is 
useful to have a concrete physical system in mind when 
considering the expected gravitational waveform. 

That pulsars are rapidly rotating neutron stars is now 
well established |22]] . Their high densities and strong 
gravitational fields allow them to withstand rotation 
rates of hundreds of times per second. Moreover, pulsar 
emission mechanisms require large magnetic fields, frozen 
into (co- rotating with) the neutron star. Indeed these 
large field strengths may produce non-axisymmetric de- 
formations of the pulsar. However, the most remarkable 
feature of pulsars is the very precise periodicity of ob- 
served pulses. 

There are more than 700 known pulsars, all at galac- 
tic distances, concentrated in the galactic plane. Based 
on the sensitivity limits of radio observations the total 
number of active pulsars in our galaxy is estimated to be 
more than 10 5 SPI- 



A. Spindown 

Pulsars lose rotational energy by electromagnetic brak- 
ing, the emission of particles and, of course, emission 
of gravitational waves [£5 26 1. Thus, the rotational 
frequency is not completely stable, but varies over a 
timescale r which is of order the age of the pulsar. Typ- 
ically, younger pulsars (with periods of tens of millisec- 
onds) have the largest spindown rates. Figure shows 
the distribution of rotational frequencies and spindown 
age, r = f/(df/dt). 
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FIG. 3. Gravitational wave frequency versus spindown 
age, t — f/(df/dt), measured in years, for 540 pulsars which 
have measured period derivative. The figure clearly shows 
a large concentration of pulsars in the mid-left of diagram. 
Most of these are isolated pulsars. The standard evolutionary 
scenario suggests that pulsars move from higher frequencies 
and shorter spindowns left and up towards this main bunch. 
In contrast, many of the millisecond pulsars lying in the upper 
right of the figure are in binary systems, and it is widely 
believed that these are pulsars which have been spun up by 
mass accretion from the companion star. 

Current observations suggest that spindown is primar- 
ily due to electromagnetic braking; however, for detection 
purposes it is necessary to construct a sufficiently general 
model of the frequency evolution to cover all possibilities. 
For observing times i bs, say, much less than r, the fre- 
quency drift is small and the rotational frequencyfj can 
be modeled as a power series of the form 



/(*) = (/o/2)(i + V/ fe * fc ) 



k 



(2.1) 



If T m in is the shortest timescale over which the frequency 
is expected to change by a factor of order unity, the co- 
efficients satisfy 



l/fcl rv r min 



(2.2) 



Clearly, for an observation time i bs <C T m i n , the first few 
terms in this series will dominate. 

Observations suggest that pulsars are born in super- 
nova explosions with very short periods (perhaps several 



3 We choose to parameterise the frequency by what will be 
the gravitational wave frequency, /o, thus introducing the ex- 
tra factor of 2 into this expression 



milliseconds), and subsequently spin down on timescales 
comparable to their age. Supernovae are observed in 
galaxies similar to our own at the rate of two or three per 
century, so we might expect r m i n ~ 40 years for pulsars in 
our galaxy. It is at this point that the distinction between 
various classes of pulsars becomes important. The known 
millisecond pulsars are old neutron stars which have have 
been spun up to periods of only a few milliseconds, possi- 
bly by episodes of mass transfer from a companion star. 
As seen from Fig. |3j, timing measurements of millisecond 
pulsars yield very long spindown timescales, T m j n > 10 7 
years. 



B. Proper Motions 

Pulsars are generally high velocity objects |2q] , as can 
be inferred by the distance they move in their lifetimes. 
Proper motions cause Doppler shifts in the observed pul- 
sar frequency. If the motion is uniform (constant veloc- 
ity), it simply induces a constant frequency shift — an 
effect which is undetectable. However, acceleration and 
higher order derivatives of the source's motion will mod- 
ulate the observed frequency. 

Studies of millisecond pulsars in globular clusters have 
shown that acceleration in the cluster field can produce 
frequency drifts which are comparable in magnitude to 
the spindown effects |27],g8| . Once again, we expect these 
effects to be well modelled by a power scries in </T cross , 
where r cross is the time it takes the pulsar to cross the 
cluster. We expect that T m ; n < r cro ss for these objects 
(since if not, the pulsar will already have escaped the 
cluster) . Thus the frequecy model adopted above should 
be sufficiently general to encompass the observational ef- 
fects of proper motions of the sources. 

A large proportion of millisecond pulsars are also in 
binary systems. Unfortunately, such pulsars participate 
in proper motions which vary over very short timescales 
(their orbital periods). The time-dependent Doppler ef- 
fect due to these motions is not mod elled accurately by 
a simple power series as in Eq. (2.1). They would re- 
quire a more elaborate model involving as many as five 
unknown orbital parameters. Including these effects in 
a coherent, all-sky search strategy would be prohibitive 
(see section |V|). In a search for gravitational waves from 
a known binary pulsar, however, it would be important 
to deal with this effect. 

Proper motions can also affect a search if the star 
moves across more than one resolution element on the 
sky during an observation. For the lengths of observation 
periods envisioned here, this is unlikely to be a problem. 
In an observation lasting a year, however, a pulsar with a 
spatial velocity of 1 x 10 3 km s _1 at a distance of 300 pc 
will move by about half an arcsecond, which is compa- 
rable to the resolution limit for our observations if the 
pulsar frequency is 1kHz. 



C. Glitches 

In addition to gradual frequency drifts due to spin- 
down, some young pulsars exhibit occasional, abrupt in- 
creases in frequency. The physical mechanism behind 
these frequency glitches is not well understood, although 
the number of observations of glitch events is grow- 
ing (£3). Given the stochastic nature of glitching, and 
the expectation that several months will elapse between 
major events, we will ignore glitching in this paper. 



III. GRAVITATIONAL WAVES FROM PULSARS 

In order to gain insight into the detection problem it is 
also important to understand the expected gravitational 
wave signal. Several mechanisms have been discussed 
in the literature which may produce non-axisymmetric 
deformations of a pulsar, and hence lead to gravitational 
wave generation ||-|l|,|l6],[||3(J . 

In general, a pulsar can radiate strongly at frequencies 
other than twice the rotation frequency. For example, a 
pulsar deformed by internal magnetic stresses, which are 
not aligned with a principal axis, can radiate at the rota- 
tion frequency and twice that frequency [pl[ . If the star 
precesses, it will radiate at three frequencies: the rota- 
tion frequency, and the rotation frequency plus and mi- 
nus the precession frequency |lq| . The important point, 
however, is that the signal at the detector is generally 
narrow band, exhibiting only slow frequency drift on ob- 
servational timescales. 

Therefore, in this section we outline the main features 
of the expected waveform and the corresponding strain 
measured at a detector for the case of crustal deforma- 
tion; other scenarios give similar results except for the 
presence of more than one spectral component. 



*xa 



L mi 



(3.4) 



is the gravitational ellipticity of the pulsar. The distance 
to the source is r, and Ijk is its moment of inertia tensor. 
The strength of potential sources is best discussed 
in terms of the characteristic amplitude h c , defined in 
Eq. (50) of 0, and simply related to h by 



(3.5) 



h c = \ hn 

\l 15 o 



For a typical IAMq neutron star, having a radius of 
10km and at a distance of lOkpc, the dimensionless am- 
plitude is 



7.7 x 10 



-25 



I zz lOkpc / /o 



10-5 10 45 gem 2 r V lkHz 



(3.6) 



The magnitude of the gravitational ellipticity, e, repre- 
sents the central uncertainty in any estimate of gravi- 
taional waves from pulsars. The tightest theoretical con- 
straint, e < 10~ 5 , is set by the maximum strain that 
the neutron star crust may support H. It has also been 
suggested that stresses induced by large magnetic fields 
might result in significant gravitational ellipticity. Re- 
cently, Bonazzola and Gourgoulhon [[L2[ have considered 
this possibility, finding discouraging results; their calcu- 
lations indicate 10~ 13 <• e <■ 10~ 9 depending on the pre- 
cise model they consider. 

In any case, an upper bound on the gravitational el- 
lipticity is e ~ 10 -5 , although typical values may be sig- 
nificantly smaller. 



B. Signal at the detector 



A. Waveform 

Adopting a simple model of a distorted pulsar as a 
tri-axial ellipsoid, rotating abo ut a principal axis with 
a frequency given by Eq. (2.1), one may compute the 
expected gravitational wave signal using the quadrupole 
formula. The two polarizations are 

h+ = h (l + cos 2 1) cos{27r/ [t + E/fe|S-]} ! (3.1) 



h x =2h cosi sm{2Trfo[t + Y,fk i k+T]} 



(3.2) 



where i is the angle between the rotation axis and the 
line of sight to the source. The dimensionless amplitude 

is 



, 8tt 2 G J Z2 / 2 

n-a - — i e 



(3.3) 



Observing the gravitational waves using an earth- 
based interferometer introduces two further difficulties 
into the detection process: Doppler modulation of the ob- 
served gravitational wave frequency, and amplitude mod- 
ulation due to the changing orientation of the detector. 

For the purpose of detection, the Doppler modula- 
tion of the observed gravitational wave frequency, due 
to motion of the detector with respect to the solar sys- 
tem barycenter, i s a l arge effect. Assuming the intrinsic 
frequency model (2.1) for the pulsar rotation, the gravi- 
tational wave frequency measured at the detector is 



/gw(t) = /o(l + ^/^ fc )( 1 



(3.7) 



where 



where v(t) is the detector velocity and h is the unit vector 
pointing to the pulsar, in some inertial frame. We gener- 
ally choose this frame to be initially comoving with the 
Earth at t = 0. The frequency measured in this frame is 



identical to that measured at the solar system barycenter 
except for an unimportant constant shift in fg. 

To understand the amplitude modulation we must in- 
troduce the Euler angles, {O, $,<!'}, which specify the 
orientation of the gravitational wave frame with respect 
to the detector frame. The dimensionless strain at the 
detector is 



h = F + (B, $, *)/i+ + F x (6, $, *)/i> 



(3.8) 



where F + and F x are the detector beam patterns given 
in Thorne M. In searching for continuous gravitational 
waves from a particular direction, the Euler angles be- 
come periodic function of sidereal time, thus resulting 
in an amplitude and phase modulation of the observed 
signal M,n2LfL9[. For observation times longer than one 
sidereal day, the amplitude modulation effectively aver- 
ages the reception over all values of right ascension, and 
over a range of declination which depends one the precise 
position of the pulsar. In particular, the effect of this pro- 
cess is to allow detection of continuous waves from any 
direction, but at the cost of reducing the measured strain 
(see Fig. 0). 



C. Parameter space 

To facilitate later discussion it is useful to parameterize 
the gravitational waveform by a vector A = (A , A) such 
that 

(\ ,X 1 ,...X 8+2 ) = (f ,n x ,n y ,f l ,...J s ). (3.9) 

Here s is the maximum number of spindown param- 
eters included in the frequency model determined by 
Eq. (£JJ) . These vectors span an s + 3 dimensional space 
on which A Q can be thought of as coordinates. (Note 



that nt 



is not an independent parame- 



ter.) In particular we denote the observed phase of the 
gravitational waveform by 



cf>(t;X) = 2irfdt / f s 



'(f) , 



(3.10) 



where f gw (t') is given by Eq. (3.7). 

Initial interferometers in LIGO should have reasonable 
sensitivity to gravitational waves with frequencies 



/ > 40Hz , 



(3.11) 



while advanced interferometers are expected to have im- 
proved sensitiviy down to 



/ > 10Hz 



(3.12) 



Moreover, theoretical constraints suggest that pulsars 
with spin periods significantly smaller than one millisec- 
ond are unlikely. This helps to constrain the highest 
frequency that one may wish to consider in an all sky 



search to be about 2kHz. According to the discussion in 
section M, the spindown parameters satisfy 



- fe < f k < T 



k 
min ' 



(3.13) 



where T m ; n is the minimum spindown age of a pulsar to 
be searched for. Finally, n x and n y are restricted by the 
relation 



X<i. 



(3.14) 



IV. DATA ANALYSIS TECHNIQUE 

Radio astronomers are familiar with searching for 
nearly periodic sources in the output of their detec- 
tors 28 32]. The technique employed by them is directly 
applicable to the problem at hand |19|j20|| . 

In the detector frame the gravitational wave signal can 
be written as 



h(t; A) = Re Ae 



(t;A) 



(4.1) 



where A — (/io+ + *^ox), ^o+ = F + (l + cos 2 i) ho and 

/iox = 2 f x (co st) hp. The orbital phase <f)(t;\) is given 

by Eqs. ( |3 . 1 0| ) and (|3.7| ). Introducing a canonical time 



t b [t;X] = 



4>{t-A) 

2vr/ 



(4.2) 



the above signal becomes monochromatic as a function 
of tjj. (The presence of the amplitude modulation com- 
plicates the following analysis without changing the con- 
clusions significantly; therefore, we treat A as constant 
in this and the next section^].) Figure [I| shows the nor- 
malised power spectr um computed from the signal as a 
function of t in Eq. ( [4.l| ) (with fk = 0), compared with 
the spectrum from the signal as a function of if,. It is 
clear that the maximum power per frequency bin is sig- 
nificantly reduced when frequency modulation is not ac- 
counted for. 



4 Amplitude modulation can be viewed as the convolution 
of the exactly periodic signal with some complicated window 
function. Thus, in reality, the power spectrum of a stretched 
signal will not be a monochromatic spike at a single frequency, 
but will be split into several discrete, narrow spikes spread 
over a bandwidth Sf ~ 10~ 4 Hz. After a preliminary detec- 
tion, the amplitude modulation spikes would provide a dis- 
criminant against false signals [nSl. 
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FIG. 4. Power spectra for two signals, each with gravita- 
tional wave frequency 5Hz, computed using approximately 10 
days worth of data; they are normalised with respect to the 
maximum power achieved if the source was directly above an 
interferometer which remained stationary during the entire 
observation. The signal was assumed to come from declina- 
tion 0° and right ascension 90°; in fact the amplitude modu- 
lation is only sensitive to changes in declination. The detector 
latitude was chosen to coincide with LIGO detector in Han- 
ford Washington. The solid line corresponds to a Doppler and 
amplitude modulated gravitational wave signal. The dashed 
line is the same signal but with the Doppler modulation re- 
moved by stretching. The (unreasonably) low frequency was 
chosen for illustrative purposes, so that both curves could 
appear on the same scale. For realistic gravitational wave fre- 
quencies (~ 500Hz) the Doppler modulated signal would be 
further reduced by roughly two orders of magnitude. 

Radio astronomers refer to this technique of introduc- 
ing a canonical time coordinate as stretching the data. 
Since interferometer output will be sampled at approxi- 
mately 16kHz, in a practical search for pulsars up to 2kHz 
gravitational wave frequency, the stretching can proba- 
bly be achieved by resampling the data stream appropri- 
ately. This method, which is called stroboscopic sampling 
by Schutz p8[ , has the benefit of keeping the computa- 
tional overhead introduced by the stretching process to 
a minimum. We will return to this issue in a later pub- 
lication. 

Now, a search of the detector output, o(t), for grav- 
itational waves from a known source is straightforward. 
One a ssum es specific parameter values £ in the wave- 
form (4.1), computes the demodulated time function 
£&[£;£] using Eq. (4.2) and stretches the detector output 
accordingly, thus 



Ob(t b [t;£]) = o(t) 



(4.3) 



If the assumed parameters £ are not too much different 
from the actual parameters A of the signal, the stretched 
data will consist of a nearly monochromatic signal. One 
then takes the Fourier transform with respect to tb, 



n^obs 



nobs 



s 2 ™^ o b (t b ) dt b . (4.4) 



Here T^ bs is length of the observation measured using t b . 
The power spectrum is then searched for excess power. 
(The threshold is set by demanding some overall statis- 
tical significance for a detection; see section W%) No- 
tice that the gravitational wave frequency, A u = /o, is 
treated somewhat differently than the other parameters; 
the Fourier transform searches over all possible values in 
a single pass. Given a sampled data set containing N 
points, the entire process, from original data through to 
the power spectrum, requires of order 3N log 2 N floating 
point operations (to first approximation). 

If all the parameters are not known accurately in ad- 
vance, it will be necessary to search over some of the 
remaining parameters A; a separate demodulation and 
FFT must be performed for each independent point in 
parameter space that one wishes to search. There are 
many possible refinements on this strategy which could 
reduce the computational cost of a search by circumvent- 
ing certain stages of the procedure described here. We 
mention some of them in section VIII, however, we focus 



attention on this baseline strategy in this paper. 

One more issue that arises in the discussion of stretch- 
ing is how it effects the noise in the detector. Throughout 
this paper we assume that the noise in the detector is a 
stationary, gaussian process; however, when we stretch 
the output data stream the noise is no longer strictly 
stationary unless it is perfectly white. Real detectors 
will have coloured noise, with correlations between points 
sampled at different times. Stretching the data modifies 
these correlations in a time dependent manner. In our 
case this is a very small effect, having a characteristic 
timescale of several hours, and besides this the noise in 
real detectors may be intrinsically non-stationary on sim- 
ilar timescales due to instrumental effects. Correcting 
pulsar searches for such non-stationarity is an important 
problem, but one that we do not address here. We sim- 
ply assume that S n (f), the power spectral density of the 
noise, can be estimated on short timescales and used in 
the conventional way for signal to noise estimates. More- 
over, the effects of stretching on noise are only a consider- 
ation when the noise is not white; since stretching affects 
the power spectrum only within bands ~ 10~ 1 Hz wide, 
the detector spectrum can usually be taken as white, un- 
less we are near a strong feature in the noise spectrum. 
The precise nature of these effects is being explored by 
Tinto §1. 



V. PARAMETER SPACE METRIC 

In general, neither the position of the pulsar nor its 
intrinsic spindown may be known in advance of detec- 
tion. Therefore, the above process, or some variant on it, 
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must be repeated for many different vectors £ until the 
entire parameter space has been explored. How finely 
must one sample these parameters in order to minimise 
the risk of missing a signal? A similar question arises 
in the context of searching for signals from coalescing 
compact binaries using matched filtering; Owen |3| has 
introduced a general framework to provide an answer in 
that case. We adapt his method to the problem at hand 
by defining a distance function on our parameter space; 
the square of distance between two points in parameter 
space is proportional to the fractional loss in signal power 
due to imprecise matching of parameters. The number 
of discrete points which must be sampled can then be de- 
termined from the proper volume of the parameter space 
with respect to this metric. 



A. Mismatch 

The one-sided power spectral density (PSD) of the de- 
tector output, stretched with parameters £ , is 



P (/)=2|5(/;£)| 2 



(5.1) 



Now, suppose a detector output consists of a signal with 
parameters A, and stationary, gaussian noise n(t) such 
that 



o(t) = h(t; A) + n(t) 



(5.2) 



Thus, the expected PSD of the detector output, once 
again stretched with parameters £, is 



E{P (f)] = 2\h(f;\,A\)\ 2 + S n (f) 



(5.3) 



where AA = £— A, and S n (f) is the one-sided power spec- 
tral density of the detector noise. (As discussed at the 
end of the previous section, we ignore the small effects of 
stretching on the noise.) The notation h(f; A, AA) indi- 
cates the Fourier transform of a signal, with parameters 
A, with respect to a time coordinate tb[t; A + AA]. We 
define the mismatch m(A, AA) to be the fractional reduc- 
tion in signal power caused by stretching the data with 
the wrong parameters, and by sampling the spectrum at 
the wrong frequency; specifically, 



m(A,AA) = 1 



|fe(/;A,AA)| s 

IM/o;A,0)| 2 



(5.4) 



Remember that A = (A = /o, A). 
In the present circumstance, it is sufficient to consider 
a complex signal 



h(t\ A) = Ae 



-2nif t h [t;\] 



(5.5) 



where the amplitude A is constant. The function tb[t; A], 
computed using Eqs. (4.2), (3.1C) and (3.7), is explicitly 
written as 



t b [t;X]= f dt'{(l + j:f k t' k )(l + v-n/c)}, (5.6) 
Jo 



Now, the Fourier transform h(f; A, AA) is 



M/;A,AA) 

where 

$[t;A,AA] 
2^ 



A 



o^obs Jo 



r obs 



dL e l ^ XA ^ 



(5.7) 



AA°4 + f (t b [t; A + AA] - t b [t; A]) (5. 



and AA° = / — fa. Here, t should be interpreted as 
a function o f tb d efin ed implicity by t b = t b \t; A + AA]. 
Using Eqs. (5.6)-( p78|) it is easy to show that m(A, AA) 
has a local minimum of zero when AA = 0; 



to(A,AA)|aa=o = 
<9aa=™.(A,AA)| aa=0 =0 



(5.9) 
(5.10) 



Thus, an expansion of the mismatch in powers of AA is 
m(A, AA) = J2 ,g a/3 (A)AA Q AA' 3 + G(AA 3 ) , (5.11) 

a, (3 



where 



9ap = \^2, dA\*d Ax em(\,A\) 

a, (3 



AA=0 



(5.12) 



In this way the mismatch defines a local distance function 
on the signal parameter space, and, for small separations 
A A, g a p is the metric of that distance function. Note that 
the metric formulation (5.11) will generally overestimate 
the mismatch for large separations, as demonstrated in 
Figure [5| 

Calculations using this formalism are considerably sim- 
plific d by partially evaluating the ri ght hand side of 
Eq. ( [5.12 ). The form of the signal ( |5.5| ) allows us to 
write 

<MA) = (d AX «$d AX0 $) - (9aa«$) <S/u/'*> . ( 5 - 13 ) 



where $ is given by Eq. ( |5.8| ), and where we use the 
notation 



1 



T'obs , 
b Jo 



(. . .)dtt 



(5.14) 



AA=0 
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FIG. 5. Fractional reduction in measured signal power 
caused by demodulating by mismatched parameters (in this 
case, an error in the assumed declination of the source). The 
solid curve is the true power ratio, the dotted is that given 
by the quadratic approximation of the metric. Note that the 
widths of the curves agree well down to 70% power reduction 
(m ~ 0.7), beyond which the metric approximation signifi- 
cantly underestimates the range of parameters permitted for 
a specified power loss. The curves are computed for a sky po- 
sition of 0° right ascension, 45° declination, and no spindown. 



B. Metric and number of patches 

Up until now, we have treated the frequency of the sig- 
nal as one of the parameters, Ao, which must be matched. 
In our search technique, stretching and Fourier trans- 
forming the data yields an entire power spectrum, au- 
tomatically sampling all possible frequencies. We would 
really like to know the number of times that this com- 
bination of procedures must be performed in a search. 
This requires knowledge of the mismatch to(A, AA) as 
a function of AA, having already maximized the power 
(i.e. minimized m) over A . The result is the mismatch 
projected onto the (s + 2)-parameter subspace: 



/' 



where 



min?7j(A, AA) 

An 



lij — 9ij 



)7«AA*AA J 



hi 



9ot90] 
goo 



(5.15) 



(5.16) 



and i — 1, . . . , s + 2. We will generally refer to /x as the 
projected mismatch. 

Technically, 7^ should be computed from g a p evalu- 
ated at the specific value of Ao at which the minimum 
projected mismatch occurred. However, since this num- 
ber is unknown in advance of detection, we evaluate "fij 
for the largest frequency in the search space. In this way 
we never underestimate the projected mismatch. 



In a search, the parameter space will be sampled at a 
lattice of points, cho sen s o that no location in the space 
has \x (given by Eq. ( 5.15 ) greater than some /z max away 
from one of the points. This is equivalent to tiling the 
parameter space with patches of maximum extent /i 
The number of points we must sample at is therefore 



1/2 



N„ 



Wdet||7ylK +2 A 



Vr 



patch 



(5.17) 



where Vp a t c h is the proper volume of a single patch, and 
s + 2 is the reduced dimensionality of the parameter space 
V (excluding Ao). 

Optimally, one should use some form of spherical clos- 
est packing to cover the space with the fewest patches. 
Our solution uses hexagonal packing in two of the dimen- 
sions and cubic packing in all the others; in this way the 
volume of a single patch is 



vpatch 



3^ /4 Mn 



4 \s + 2 



(s+2)/2 



(5.18) 



VI. DEPTH OF AN ALL SKY SEARCH 

We are finally in a position to estimate the depth of a 
search for periodic sources using LIGO. The detector par- 
ticipates in two principal motions which cause significant 
Doppler modulations of the observed signal: daily rota- 
tion, and revolution of the Earth about the Sun. The 
latter is actually a complex superposition of an ellipti- 
cal Keplerian orbit with a smaller orbit about the earth- 
moon barycenter, and is further perturbed by interac- 
tions with other planets. For now, however, we use a 
simplified model which treats both rotation and revolu- 
tion as circular motions about separate axes inclined at 
an angle e = 23°27' to each other. Although a simplifica- 
tion, this does remove any spurious symmetries from the 
model; thus, an actual search using the precise ephemeris 
of the earth in its demodulations should give comparable 
results. In this model, then, we write the velocity of the 
detector in a frame which is inertial to the solar system 
barycenter but initially comoving with the earth: 



- (QRd sin Qt — QaRa sin Q,a€)X 
-(QRd cos fit — VLaRa cos e[cos CI At 
-CIaRa sin e[cos ft At — l])z 



m (6.1) 



where R c i = 6.371 x 10 8 (cosZ)cm, I is the latitude of 
the detector, and Ra = 1.496 x 10 13 cm is the distance 
from the earth to the sun. The angular velocities are 
fl = 2tt/(86 400s) and fi A = 2vr/(3. 155674 x 10 7 s). Our 
coordinate system measures x towards the vernal equinox 
and z towards the north celestial pole, and we arbitrarily 
choose to measure time starting at noon on the vernal 
equinox. 
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The number of spindown parameters fk which must be 
included to account for all intrinsic frequency drift de- 
pends to a large extent on the type of pulsar one wishes 
to search for. We determined this number on a case by 
case basis, including all parameters which lead to a signif- 
icant increase in the number of parameter space patches. 
Equivalently, the following geometric picture suggests a 
simple criterion for deciding when there is one spindown 
parameter too many included in the signal parametriza- 
tion. Let X L be the last, 'questionable' spindown parame- 
ter f s (so L = s + 2). With respect to the natural metric 
7y on parameter space, the unit-normal to surfaces of 
constant X L is just Y L l{l LL ) 1 ^ 2 1 where 7 y is the inverse 
of jij . The spindown parameter X L is unnecessary if the 
proper thickness of the parameter space in this normal 
direction is less than half the proper grid spacing; that is, 
if 2r min +2 /(7 LL ) 1/2 < vt/i' In practice, one has in- 
cluded more spindown parameters than necessary if and 
only if 7 LL > ^r-^+V/w 



A. Patch number versus observation time 

It is extremely difficult to obtain a closed form expres- 
sion for the metric, let alone its determinant. Therefore, 
we present results for two concrete scenarios which sug- 
gest themselves based on the discussion in section |J: (i) 
hypothetical sources with / < 1000Hz, and spindown 
ages greater than r = 40Yrs; incidentally, this also in- 
cludes the majority of known, millisecond pulsars; and 
(ii) slower sources (/o < 200Hz) having spindown ages in 
excess of r = lOOOYrs. The number of parameter space 
points which must be searched is plotted as a function of 
total observation time in Fig. |6J. The numbers are nor- 
malised by a maximum projected mismatch ^ max = 0.3. 

In considering an optimal choice of observation time, it 
is useful to construct an empirical fit to -/V p (£ bs,Mmax)- 
Noti ce fi rst that all the parameters AA in $, given by 
Eq. ( |5.8| ), appear multiplied by the gravitational wave 
frequency /o; thus, N p oc (/o) s+2 - Furthermore, provided 
the determinant of the metric is only weakly dependent 
on the values of the fk one may also extract a factor of 
r -s(s+2)/2. our investigations suggest the validity of this 
approach. In this way we arrive at the expression 



and T is the observation time measured in days. These 
formulae are normalised using only the data correspond- 
ing to Fig ||(a), and subsequently compared with com- 
puted values for several frequencies and spindown ages 
t. The analytic fit is in good agreement with the com- 
puted results for a variety of parameters; however, the 
fits generally break down for observation times less than 
one day. We stress that more spindown parameters may 
become important for observation times longer than 30 
days. 

Schutz |18[ has previously estimated the number of 
points which must be searched in the abscence of spin- 
down corrections; he argued that this number scaled as 
T 4 for observation times longer than about a day. The 
difference b etwe en his previous estimate and the expres- 
sion in Eq. (3.4), which shows that the number of points 
increases as T , derives from an asymmetery between 
declination and right ascension which was not accounted 
for in his argument. 

The benefit of the metric formulation is that it ac- 
counts for the significant correlations which exist between 
the intrinsic spindown and the earth-motion-induced 
Doppler modulations by using points which lie on the 
principal axes of the ellipsoids described by Eq. (5.15) 



Replacing the invariant volume integral in Eq. (5.17)by 



IL7« d a+2 X 



(6.7) 



gives the number of points required for a search if, in- 
stead, one chooses them to lie on the {n x , n y , /i, fe, ...} 
coordinate grid. Figure ffl shows the total number of 
points computed using this method compared to the re- 
sults obtained using the invariant volume integral. For 
sufficiently long integration times the difference can be 
several orders of magnitude. 



N p ~ max [Af s F s (t)} 



(6.2) 



where 



( h 



s+2 



/40Yrs\ s(s+1)/2 / 3 ^ (s+2)/2 



fin 



V lkHz 7 V t J 

Fo(i bs) = 6.9 x 10 3 T 2 + 3.0 T 5 

1.9 x 10 8 T 8 + 5.0x 10 4 T n 



Fl(tobs) 

-FM^obs) 



4.7 + T 6 
2.2 x 10 7 T 14 
56.0 + T 9 ' 



(6.3) 
(6.4) 
(6.5) 

(6.6) 
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FIG. 6. Number of independent points in parameter space 
as a function of total observation time, using a maximum pro- 
jected mismatch /i max = 0.3. The parameter ranges chosen 
were: (a) maximum gravitational wave frequency 1000Hz, 
minimum spindown age r m m = 40Yrs (hypothetical young 
pulsars); (b) maximum gravitational wave frequency 200Hz, 
minimum spindown age Tmin = 10 3 Yrs (observed, slow pul- 
sars). The short-dashed curve represents the total number of 
patches ignoring all /*.. The long-dashed curve is the num- 
ber of patches including only /i in the search. The dotted 
line is the number of patches including both /i and /a . Also 
shown is the empirical fit given in the text; it was normalised 
by the results in shown in (a). In some regimes, searching 
over an additional spindown parameter would seem to reduce 
the number of patches; however, this actually only indicates 
regions where the parameter space extends less than one full 
patch width in the additional dimension. In such regimes one 
must properly discard the extra parameter from the search, 
forcing one to choose always the higher of the curves. 
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FIG. 7. The total number of parameter-space points 
needed to search for pulsars having gravitational wave fre- 
quency up to 1kHz, and spindown age greater than r = 40Yrs. 
The solid line is the number computed using the metric and 
properly accounting for correlations between various terms in 
the frequency evolution. The dotted line is the same num- 
ber computed directly by assuming the points must lie on the 
grid of coordinates used to parameterise the signal. The ben- 
efits of using the metric to optimally place the points to be 
searched in parameter space is clear. 



B. Computational Requirements 

The number of real samples of the interferometer out- 
put for an observation lasting i b s seconds, and sampled 
at a frequency 2/ max , where / max is the maximum grav- 
itational wave frequency being searched for, is 



ti -^/max^obs 



(6.8) 



For each A that is used to stretch the detector output, 
a search then requires an FFT, calculation of the power, 
and some thresholding test for excess power. Assum- 
ing that the stretching and thresholding require negligi- 
ble computations compared to performing the FFT and 
computing the power, the total number of floating point 
operations for a search is 

A^op = 6/ max i obs Ayiog 2 (2/ max i obs ) + 1/2] , (6.9) 



where N p is given by (6.2)-(6.6). The additive 1/2 in- 
side the square brackets accounts for the three floating 
point operations per frequency bin which is required to 
compute the power from the Fourier transform. 

A guideline for a feasible, long-term, search strategy is 
that data reduction should proceed at a rate comparable 
to data acquisition. Thus, the total computing power 
required for data reduction, in floating point operations 
per second (flops), is 
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N 

P = T^- = 6/ max iVp(i o bs,Mmax)[l0g 2 (2/ m a X £obs) + 1/2] . 



t-obs 



(6.10) 



For a prescribed maximum projected mismatch p. max , 
and maximum available computing power P max this ex- 
pression determines the maximum allowed coherent inte- 
gration time. Alternatively, given the computing power 
available for data reduction, P max , it provides an implicit 
relation between /z max and the integration time. 

The idea now is to choose /imax and t b s so that we 
maximize the sensitivity of the search. In order to do this 
we must first obtain a threshold, above which we consider 
excess power to indicate t he p resence of a signal. 

As discussed in section [TV], we assume that the noise 
in the detector is a stationary, gaussian random process 
with zero mean and PSD S n (f). In the abscence of a 
signal, the power P (f) = 2|n(/)| 2 is exponentially dis- 
tributed with probability density function 



-Po(J)/S n tf) 



(6.11) 



We assume that there is independent noise in each of 
/max^obs frequency bins for a given demodulated power 
spectrum. In general the noise spectra obtained from 
neighbouring parameter space points will not be statis- 
tically independent; however, one may expect that the 
correlations will be small when the mismatch between 
the points approaches unity. Therefore we approximate 
the number of statistically independent noise spectra in 
our search to be N p (t ba , Mmax = 0.3). In order that a 
detection have overall statistical significance a, we must 
set our detection threshold so there is less than 1 — a 
probability of any noise event exceeding that threshold. 
For a detection to occur the power in the demodulated 
detector output must satisfy 



> 



Pr 



PoU) 

Sn{f) ' S n (f) 



hi 



Jmax^obs^ *p\^obs j /'max — u.oj 



.12) 



where P (f) was defined in Eq. (IO), and p c is the thresh- 
old power. 

In other words, if the power at a given frequency ex- 
ceeds p c we can infer that a signal is present; the expected 
power in the signal is then p c — S n . Thus, the minimum 
characteristic amplitude we can expect to detect is 



tth 



( Pc /S n - l)S n (f) 



<Fl(Q,$,*)>(l-<n>)t ohB ' 



(6.13) 



where < F?(0, $, ^>) > is the square of the detector re- 
sponse averaged over all possible source positions and 
wave polarizations. < p > is the expected mismatch for 
a source whose signal parameters A lie within a given 
patch, assuming that all parameter values in that patch 



are equally likely. We note that the characteristic de- 
tector sensitivities h s / yr in Fig. |l] are obtained from this 
expression by setting £ b s = 10 7 seconds, < p >= 0, and 
/max ^obs -Mp = 1 in the expression for p c ; this agrees with 
Eq. Q. 

The optimal search strategy is to choose those values 
of i bs and /z max which, for some specified computational 
power P max and detection confidence a, maximize our 
sensitivity O which is defined by 



0(*obs,Mmax) = T— OC 

rith 



'\l_s±2 Ij, 

l ± s -|-4 MmaxJ - 1 



Pc/S n 



1 



(6.14) 



where p c /S n is given by Eq. (3.2). Assuming an overall 
statistical significance of a = 0.99, we have computed 
the optimal observation time i b s and optimal maximum 
mismatch p maDC , as functions of computing power, for the 
two searches considered in the previous subsection. The 
results are shown in Fig. pi 



VII. COMPUTATIONAL REQUIREMENTS FOR 
A DIRECTED SEARCH 



In sections M and VI we examined the computational 
requirements of an all-sky pulsar search. In this sec- 
tion we examine the computational requirements for a 
directed pulsar search, by which we mean a search where 
the position is known but the pulsar frequency and spin- 
down parameters are unknown. Obvious targets in this 
category are SN1987A, nearby supernova remnants that 
do not contain known radio pulsars, and the center of 
our galaxy. Such searches will clearly be among the first 
performed once the new generation of gravitational wave 
detectors begin to come on line. 

Our treatment of directed pulsar searches closely par- 
allels that of of the all-sky search, so we can be brief. 
Since the source position (n 



xi "y, 



is known, we can sim- 
ply remove the Earth's motion from the data. Below we 
imagine that the signal has already been transformed to 
the solar system barycenter. Then the unknown param- 
eters describing the pulsar waveform are 



(A°,A 1 ,...A s ) = (/ ,/ 1 



./.) 



(7.1) 



where the fi are the same as defined in Eq. (pTTJ) and 
s is just the number of spindown parameters included 
in the frequency model. We a gain calcu late the met- 
and 7y using Eqs. ( 5.13 ) and ( 5.1 6| ) respectively, 



and then calculate N p using ( 5.17 ) (except the inte- 



gral is now over s-dimensional parameter space). As- 
suming hexagonal packing in two dimensions and cu- 
bic packing in the others, the size of each patch is 
F pa tch = (3\/3/4)(4^ max /s)' !! / 2 . (Except for s = 1, where 
^p a tch = 2/XmaxO We arrive at the expression 



N p ~ max [Af s G s (t)] , 

s6{l,2...} 



(7.2) 
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where 



K 



( h V/40Yrs\ 
^lkHz/ V t J 
G 1 (U s ) = 1.5xl0 3 T 2 
G 2 (U s )=6.97xl0 1 T 5 
G 3 (t oha ) = 2.89 x 1(T 4 T 9 , 



«(s+l)/2 



0.3 



8/2 



(7.3) 

(7.4) 
(7.5) 
(7.6) 



where T is the observation tim e m ea sure d in days. Com- 
paring these results with Eqs. ( |6.2j )-(3.6), we see that for 
our fiducial parameter values (/o = 1kHz, r min = 40Yrs, 
A*max = 0.3) and observation times T of order a week, 
N p is ~ 1 5 times larger for an all-sky search than for a 
directed search. Another way of putting this is: after us- 
ing one's freedom to adjust the frequency and spin-down 
parameters in optimizing the fit, only ~ 10° distinguish- 
able patches on the sky remain. Equivalently, a single 
directed search can cover an area of ~ 10~ 4 steradians. 
Thus ~ 1000 week-long, directed searches would be suf- 
ficient to cover the galactic center region. 

We can calculate the optimal /i max and £ bs as a func- 
tion of computing power for a directed search in the same 
way as we did for th e all -sky directed search. (Except 
the factor 2±| in Eq. 5.14 becomes j^ for the directed- 
search case.) The results are shown in Fig. [| for our 
two fiducial types of pulsar. We see that knowing the 
source position in advance increases £ bs by only a factor 
of ~ 10, for 1 Tflops computing power. The resulting 
gains in sensitivity can be seen in Fig. @. 
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FIG. 8. The optimum observation time (thick solid line), 
and maximal projected mismatch (dotted line) as functions of 
available computational power. Both graphs assume a thresh- 
old which gives an overall statistical significance of 99% to 
any detection (although the results should be insensitive to 
the precise value). Each of the graphs corresponds to: (a) 
the situation encountered when searching for periodic sources 
having gravitational wave frequencies up to 1000Hz, with min- 
imum spindown ages r m i n = 40Yrs. (b) The equivalent results 
for gravitational wave frequencies up to 200Hz, with minimum 
spindown ages r m i n = 10 Yrs. The transition region seen in 
figure (a) is due to the fact that a longer integration time 
would require searching over an additional spindown param- 
eter, as seen in Fig. H. In this region it is more efficient, as 
one adds computational power, to lower mismatch thresholds, 
rather than searching over the additional parameter. 
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FIG. 9. The optimum observation time (thick solid line), 
and maximal projected mismatch (dotted line) as functions 
of available computational power for directed searches. Both 
graphs assume a threshold which gives an overall statistical 
significance of 99% to any detection (although the results are 
insensitive to the precise value). Each of the graphs corre- 
sponds to: (a) the situation encountered when searching for 
periodic sources having gravitational wave frequencies up to 
1000Hz, with minimum spindown ages r m in = 40Yrs. (b) 
The equivalent results for gravitational wave frequencies up 
to 200Hz, with minimum spindown ages r m i n = 10 3 Yrs. The 
transition regions, where the optimum observation time does 
not increase, are due to the fact that a longer integration time 
would require searching over an additional spindown param- 
eter. 



VIII. FUTURE DIRECTIONS 

Searching for unknown sources of continuous gravita- 
tional waves using LIGO, or other interferometers, will be 
an immense computational task. In this paper we have 
presented our current understanding of the problem. By 



applying techniques from differential geometry we have 
estimated the number of independent points in the pa- 
rameter space which must be considered in all-sky and di- 
rected searches for sources which spin down on timescales 
short enough to produce observable effects; these num- 
bers were used to compute the maximum achievable sen- 
sitivity for a coherent search (see Fig. 0). Furthermore, 
the metric formulation can be used to optimally place 
the parameter space points which must be sampled in a 
search. 

Our analysis takes no account of bottlenecks in the 
analysis process due to data input/output and inter- 
processor communication. These are important issues 
which may impose further constraints on the maximum 
observation time; however, it seems premature to address 
such problems until we know the hardware that will be 
used to conduct searches for continuous waves. 

Unfortunately, Fig. @ shows that it will be impossible 
to search, in one step, 10 7 seconds worth of data over all 
sky positions. However it is also unnecessary. We fore- 
see implementing a hierarchical search strategy, in which 
a long data stream is searched in two (or more) stages, 
trading off sensitivity in the first stage for reduced com- 
putational requirements. Having determined a number 
of potential signals in the first stage — presumably at a 
threshold level which allows many false alarms due to 
random noise — these candidate events would be followed 
up in the second stage, using longer integration times. 
The longer integration times would be possible because 
the search would only have to be performed over much 
smaller regions of the parameter space, in the neigh- 
bourhoods of the candidate signal parameters. In this 
way, one can achieve a greater sensitivity than a coher- 
ent search using the same computational resources. 

Clearly one can imagine many different implementa- 
tions of this rough strategy, and we have not yet deter- 
mined the optimal one. Nevertheless, we have consid- 
ered the simple example where the data is searched in 
two stages. Candidate signals from an all-sky search of a 
short stretch of data [T^ seconds long] are followed up 
using longer Fourier transforms to achieve greater sensi- 
tivity. One can estimate T^ using Fig. || and an assump- 
tion that roughly half of the total computing budget is 
used on the first stage; this turns out to be a valid as- 
sumption. A simple argument along these lines goes as 
follows. Consider a search for 'young, fast' pulsars that 
begins by coherently analyzing stretches of data that are 
all ~ 1 day long (possible with ~ 4x 10 12 flops, by Fig. ||). 
Imagine that in the second stage of the search one fol- 
lows up all templates such that Po(f, A) > 4.6S n (f), by 
seeing whether templates with roughly the same parame- 
ter values are exceeding this threshhold every day. (Here 
Po(f, A) is the power of the stretched data at frequency /, 
for stretch A. This threshhold implies that one is follow- 
ing up only one out of every hundred templates.) It seems 
likely that this second stage will not be more computa- 
tionally intensive than the first. To exceed this thrcsh- 



17 



old, a pulsar must have h c > 12/i 3 / yr . This is factor of 
roughly 3 better than if one restricted oneself to coherent 
searches considered above, but is a factor of 3 worse than 
the sensitivity one could achieve with unlimited comput- 
ing power. 

A refinement of this strategy would be one in which 
the first pass consists of several incoherently-added power 
spectra. That is, one slices the data into N sequential 
subsets, performs a full search (as described in this paper) 
for each subset, and adds up the power spectra of the re- 
sulting searches for each of the parameter sets. This tech- 
nique has been used to good effect by radio astronomers 
searching for pulsars p8[ . Since the addition of power 
spectra is incoherent, there is a loss of signal-to-noise ra- 
tio in the final summed power spectrum of 1/yN in re- 
lation to a full coherent search over the whole timescale. 
However, the computational savings involved allow one 
to search stretches of data which are much longer over- 
all. For some optimal choice of N, this will result in 
higher sensitivities when one follows up candidate de- 
tections using coherent searches. D. Nicholson (private 
communication) has estimated that a lTflops computer 
could perform such a search of 10 7 s of data, over all sky 
positions but ignoring pulsar spindowns. A subsequent 
paper will present a concrete analysis of this and other 
hierarchical scenarios |34[ |. 
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